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Errors Introduced by Finite Space and Time Increments 
in Dynamic Response Computation 

Samuel Levy and Wilhelmina D. Kroll 

An investigation is made of the accuracy and stability of numerical integration methods 
when applied to the computation of the dynamic response of structures to impact loads. 
The effect of finite time increments is studied both by obtaining analytical solutions for a 
single-degree-of-freedom system and by carrying out numerical integrations for many- 
degree-of -freedom systems; the effect of finite space increments is studied by replacing a 
continuous beam by a discrete number of elastically connected point masses. It is found that : 
(1) Of the methods investigated, only Houbolt's is stable when the time increments are largo 
compared with the natural periods of the system. Errors are introduced by Houbolt's 
method, in this case, which result in the damping out of the responses in the higher modes 
of vibration. All of the methods give good results when the time increment is less than about 
1/30 of the period in the highest frequency mode. (2) The distributed mass of a beam can 
be considered to be concentrated at relatively few mass points for computational purposes; 
using a five mass idealization, the bending moment at the center of a uniform beam is deter- 
mined with good accuracy. 



1. Introduction 

With larger aircraft, higher landing speeds, and the 
necessity of flying in bad weather, the transient 
vibrations caused by severe gusts, landing impacts, 
and similar shock loads are becoming increasingly 
important in the stress analysis of airplanes. The 
usual method of computing these transients is to 
superpose the response in a small number of the im- 
portant modes [1, 2]. 1 Such computations are 
lengthy and in some cases give results that are of 
questionable accuracy [3]. 

Houbolt [4] presents a numerical integration 
method in which the derivatives in the equations of 
motion are replaced by finite differences to permit a 
step-by-step calculation of the dynamic response of 
an elastic aircraft entering a gust. Houbolt's 
method is adaptable to the problem of the dynamic 
response of an airplane to landing impact and was 
used [5] to determine the deflections of the wing for 
an unsymmetrical two-point landing of a model air- 
plane. From these deflections, the bending mo- 
ments on each wing at stations 17.5 inches from the 
wing root were computed. The computed results 
were compared with experimental results. The 
agreement was good, indicating that the difference 
equation approach holds promise as a means of 
determining the dynamic response of airplanes to 
shock loads. 

The purpose of the present paper is to determine 
the effect of the use of finite time increments and of 
the replacement of the continuous structure by a dis- 
crete number of elastically connected point masses 
on the accuracy and convergence of numerical inte- 
gration methods. The errors due to using too 
coarse a time increment to adequately describe the 
fine detail in the force-time history and due to 
approximating the initial conditions by various finite 
difference approximations are not considered in this 

1 Figures in brackets indicate the literature references at the end of this paper. 



report. A general study of errors in numerical inte- 
gration procedures, using techniques similar to those 
in this report, is given in references [6] and [7]. 

2. Error Due to Finite Time Increment 

In order to determine the errors introduced by 
finite time increments in numerical integration meth- 
ods, a study is made of a single-degree-of-freedom 
system having simple loading conditions. The mo- 
tion of more complicated systems can be considered 
as being made up of the motion in several normal 
modes, each mode acting as a single-degree-of-free- 
dom system. The results of this study can be used, 
therefore, in judging the adequacy of the various 
numerical integration methods for complicated sys- 
tems as well as for simple systems. The results are 
determined in analytical form, using the calculus of 
finite differences, for convenience in judging the 
accuracy of the methods and the peculiar nature 
of the errors introduced by using them. 

2.1. Basic Problem Considered 

The basic problem considered in this section is 
shown in figure 1. The mass, if initially disturbed, 
should vibrate without damping or amplitude build- 
up at a natural frequency u=^k/m (rad/sec). For 
such a system, the equation of motion is 



m -rp-\-kx = 0. 



(i) 



In applying numerical integration methods to the 
solution of this equation, the following notation is 
used: 

A^ — time increment between successive steps in 
the numerical integration process. 
/i = number of steps taken from £=0 to t=nAt; 
as subscript, indicates value when t = nAt. 

x n = displacement when t=nAt. 
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Figure 1. Basic problem considered in comparing numerical 
integration methods. 



2.2. Method Replacing Second Derivative by 
Equivalent Central Differences 



Substituting into eq (1) 
1 



\dt 2 



).■ 



w 



\ x n+\ 2j , n + J' w _i) 



gives 



m 



/ A A2^» + 1 2£n\ x n- \)\lCXn — 0, 



which reduces to 



-+FM 



x n +x n -i=0. 



(2) 



(3) 



(4) 



To solve eq (4) by the ealcxilus of finite differences, 
we make the substitution, 



x n =Ap n , 



(5) 



where A is the arbitrary constant to be determined 
from initial conditions, and fi is the number to be so 
chosen that eq (4) is satisfied. Substituting from 
eq (5) into eq (4) gives 



^/3 n+1+ r^^ 2 ~ 2 l A/3 n + At3 n - 



-0. 



Dividing through by Afi n \ eq (6) reduces to 



3 2 + 



PM 



1+1 = 0. 



(6) 



(7) 



Several cases are of interest: 

Case 2.2.1 : 0<A£<2 A /m/&; we will first consider 
At=^2^m/k; that is, At=0.225(2ir)-sm/k, where 



(27r)^m/k is the natural period of the system. In this 
case, eq (7) reduces to 



£ 2 +l=0, 



from which, 

t3=±^-A = 
Substituting this value of into eq (5), 



■±i=« ±2 



(8) 
(9) 



x n = Ae 2 +Be ~ 2 =A' sin 2 +B' cos 2 , (10) 

where A, B, A' , B' are arbitrary constants deter- 
mined from the initial conditions. Since 

n=t/(At)=t^Jk/m/^2 f (11) 

eq (10) can be written 

x=A' sin l.llt^Tcfm + B' cos l.llUjk/m. (12) 

In this case, the effect of the numerical integration 
method is to increase the effective natural frequency 
from v^/m to l.ll\'k/m without introducing damping 
or buildup of the response. 

As the value of the finite time increment At is 
varied from to 2-y/m/k, the factor 1.11 in eq (12) 
varies from 1 to x/2, but otherwise the form of eq 
(12) does not change. 

Case 2.2.2: 2^/m/k<^At; we will first consider 
At=3^m/k. Although such a large value of At 
would not ordinarily be used on a single-degree-of- 
freedom system, it is usually impossible to avoid such 
a large value in a many-degree-of-freedom system 

having high frequency modes. With At=3^m/k, 
eq (7) becomes 

£ 2 +7/3+l = (13) 

giving 

/3=-.1459 = -<r 1 - 925 ; 0= -6.8541 = -e 1 - 926 . (14) 

Substituting these values of (3 into eq (5), 

x n =A(—l) n e' 1MH +B(—l) H e lMn 

=A' cos mr sinh 1. 925ft +5' cos nw cosh 1.925ft. 

(15) 

Substituting n=tjAt={t-\k/m)/ , 3 into eq (15) gives 

x=[A' sinh 0.642*V*/ra + 5' cosh 0.642/y£/m] 
cos 1 .047 t\kjm. (16) 

In this case, the primary effect of the numerical 
integration method is to introduce hyperbolic 
functions of time in the answer. Because these 
functions increase indefinitely with time, the solution 
is unstable. 
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For all values of the time increment AC>2^mjk, 
unstable solutions of the type shown by eq (16) re- 
sult. The larger the value of At, the greater will be 
the instability. It is this fact that makes the method 
considered in this section (2.2) inappropriate for use 
on many-degree-of-freedom systems where it is diffi- 
cult to avoid large values of At relative to the higher 
frequency modes. 

2.3. Houbolt's Method 
If, in eq (1), Houbolt's substitution 

/ d 2 x\ 1 

\jf) =T£jx2(2Xn--5X n - 1 + 4x n _.2--Xn-3) (17) 

is made, we obtain 

Tfl 

(Aft 1 ( 2Xn ~ 5 - r »-i + 4- r »-2 — .r„_ 3 ) + £.r„ = 0. (18) 
After transposing terms, eq (18) can be written as 



2 + 



m 



x n — 5x n -i+4x n - 2 —Xn-3=0. (19) 



Making use of the substitution given in eq (5), eq 
(19) becomes 

L M J (20) 

Dividing through by Afi n ~ 3 gives 



D + r> 



5j8 2 + 40-l = O. (21) 



Several cases are of interest. 



Case 2.3.1: We will again consider M= y2 -y/m/Jc, 
as was done in case 2.2.1. Equation (21) then 
reduces to 

40 3 -5/3 2 + 4/3-l = O, (22) 

giving the roots, 

£=0.3710 = £-° 991 \ 

£=0.439 + 0. 694i = £ ( -°- 197+1 - 007i) > (23) 

/3 = 0.439-0.694i = f ( -°- 197 - 1 - 007 ^. ) 
Substitution of these values of f3 into eq (5) gives 

x -_^4^-0.991n_|_2Jg(-0.197« + 1.007in) _!_ (J e (-0.197n -1.007*n) 
= ^-0.991n + g-0.197»( S / gm 1.007ft + C" COS 1.007/?) 

_ (24) 

Using n = tlM={t^kjm)l^2 gives 

( B' sin 0.712* J— 



x=Ae-°™ 1 * <Jk/m-\- e -o.iMtjk/m [ B f sin 0.712J +1- 

117 \ 
C cos 0.712^ 



In this case, the effects of the numerical integration 
method used are (1) to introduce a decaying ex- 
ponential, the first term; (2) to lower the natural 
frequency, 0.712 instead of 1; and (3) to introduce a 
damping term on the response, e-° 139 <V*7^. 

Case 2.3.2. The errors in the case where A is 
small can be judged best by a numerical example. 
With At=^0.l^m/k, a solution similar to that for 
case 2.3.1 gives 



x = A e - 2 - 3 2t V*/™ -I- e -00095* yV 
C cos 0.961* 



(lT rin 0.961* ^+ 

V£) (26) 



In this case, the errors are a small drop in natural 
frequency, 0.961 in place of 1; and the presence of 
decaying exponentials. In the limit when Af->0, 
the 'coefficient 2.32^oo; the coefficient 0.0095-K); 
and the coefficient 0.961-^1. The solution ap- 
proaches the exact solution as At approaches zero. 

Case 2.3.3. To judge the errors for large values of 
A*, we will consider the numerical case where At= 
V98ym/&. In this case, the solution is given by 



C cos 0.183/ J^V 



//r sin 0.183*-%/ 



k 



+ 



(27) 



The errors are a large drop in natural frequency, 
0.183 in place of 1, and the presence of decaying ex- 
ponentials. In the limit when At-^co, the coeffi- 
cients 0.1805, 0.1424, and 0.183 all approach zero. 
It is thus seen that large errors may result for large 
values of At, but that they will always be of a stable 
type; rounding-off errors will not build up. 

3. Error Due to Replacing Continuous Struc- 
ture by a Discrete Number of Elastically 
Connected Masses 

An airplane wing is, essentially, a beam and, like 
a beam, has an infinite number of degrees of freedom. 
Since an analysis of such a structure without simplifi- 
cation be difficult, if not impossible, it is customary 
to replace the distributed mass of the airplane by 
lumped masses connected by massless springs chosen, 
appendix II of reference [8], to approximately repre- 
sent the large-scale motions of the airplane. It is 
also customary to approximate the distributed mass 
by generalized masses [1], corresponding to the lower 
modes of vibration. Such lumped or generalized 
masses are never equivalent in their behavior to the 
distributed mass of the airplane, since they can only 
represent a finite number of degrees of freedom. 

A general investigation of errors resulting from the 
use of a finite number of masses is beyond the scope 
of this report; however, since the Houbolt method 
has advantages insofar as numerical integration is 
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m> ej_ 



Figure 2. Beam used in computations. 



t 

Figure 3. Time history of force applied to beam. 

concerned, a study is made of the error resulting 
when this method of numerical integration is used 
and a finite number of lumped masses replace the 
distributed mass. This is done by comparing the 
computed responses with 3, 5, and 7 lumped masses 
for an undamped, uniform, free-free beam of mass 
ra, length 21 and bending stiffness, EI, figure 2, using 
Houbolt's method as modified in reference [5]. It is 
assumed that the beam is subjected to suddenly ap- 
plied transverse loads of the type shown in figure 3. 
As more complicated time histories of loading can be 
considered to be made up of stepwise changes of the 
type shown in figure 3, the results will be significant 
for any type of transverse loading. All displace- 
ments are taken in a vertical direction. Translation 
and rotation of the beam are considered in addition 
to bending. 

3. 1. Uniform Beam With Load Applied at Center 

A uniform beam is considered to be suddenly sub- 
jected to a constant normal force P at its center. 
The beam is idealized by considering the mass of 
the actual beam, figure 4 (a), to be concentrated at 
3-, 5-, and 7-mass points, figures 4 (b), 4 (c), and 
4 (d), respectively. 

To evaluate the combined errors due to the use 
of finite numbers of mass points and finite time in- 
crements, three time increments are used in the 
numerical analysis for each mass distribution. 
These time increments are 

At=^mFJl2EI, ^mP/l2EIf and ^mP/12EI- 



m, £f 



(a) Actual beam 



m/4 m/2 m/4 



(b) Three mass system 
m/8 m/4 m/4 m/4 m/8 

oa^a-Q^vvmO-^vvhOaaa-o 

(c) Five mass system 



m/12 m/s m/6 m/6 m/6 m/6 m/12 

0-na/\-0 _vv v-0-^ /V vO~ vv ^^ 

(d) Seven mass system 

Figure 4. Distributions of mass considered for beam subjected 
to impact at center. 



The computations are carried out as described in 
reference [5] and as shown for a typical case in the 
appendix. In this procedure, the initial conditions 
are taken into account by considering the displace- 
ments zero for three time increments prior to the 
initial instant of time. The flexibility is taken into 
account by the use of influence coefficients. Matrix 
algebra is used to simultaneously compute the dis- 
placements at all the mass points at a particular in- 
crement of time from the given forces acting at that 
time and from the known displacements at the pre- 
ceding three increments in time. The bending 
moment at the center is then computed from the 
inverse of the influence coefficient matrix, the dis- 
placements at the various mass points, and the lever 
arms from these mass points to the center of the 
beam. The bending moment, M, at the center of 
the beam, positive when concave away from the side 
of the beam being struck, is plotted on a dimension- 
less basis in figures 5, 6, and 7 for the 3-, 5-, and 7- 
mass approximations, respectively. It is noted that, 
as the time interval used in the step-by-step nu- 
merical integration is shortened, the maximum abso- 
lute value of the bending moment is increased. As 
is to be expected from the analysis of a single de- 
gree-of- freedom system, the errors of the numerical 
method cause a damping out of the oscillatory re- 
sponse and, as the time increment is decreased, a 
decrease in the period of the oscillatory response. 

For each of figures 5 to 7, the mass distribution is 
kept constant, and the increment of time used in 
the numerical integration is varied to show the effect 
of time increment. In each of figures 8 to 10, the 
time increment is kept constant, and the responses 
for the different mass distributions are plotted to 
bring out more clearly the effect of increasing the 
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Figure 5. Bending moment ratio at the center of the beam 
caused by impact at center, three-mass idealization. 
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Figure 6. Bending moment ratio at the center of the beam 
caused by impact at center, Jive-mass idealization. 
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Figure 7. Bending moment ratio at the center of the beam 
caused by impact at center, seven-mass idealization. 
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Figure 8. Bending moment ratios for different mass 
idealizations with impact at center mass, time increment 
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Figure 9. Bending moment ratios for different mass 
idealizations with impact at center mass, time increment 
At=l/3^mP/l2Er. 
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Figure 10. Bending moment ratios for different mass 
idealizations with impact at center mass, time increment 
At=lj^miyi2EL 

number of masses. The maximum absolute value 
of the bending moment at the center differs by not 
more than 4 percent for the 5- and 7-mass systems, 
regardless of time increment. There is some differ- 
ence in the period of the response for the shorter 
time increments for the 5- and 7-mass systems. 

3.2. Beam With Impact at One Tip 

Because it was believed that an impact at a tip 
mass might result in a more severe loading condition, 
an investigation is made for the case of a load P 
applied at the tip of the beam shown in figure 4 (a), 
using the idealizations shown in figures 4 (b) and 
4 (c). The bending moments at the center for the 
3- and 5-mass idealizations are shown in figures 1 1 
and 12, respectively. It is noted again, as in case 
3.1, that as the time interval used in the step-by- 
step numerical integration is shortened, the maximum 
absolute value of the bending moment at the center 
is increased. As is to be expected from the analysis 
of a single-degree-of-freedom system, the errors of 
the numerical method cause a damping out of the 
oscillations, which is particularly evident for the 
longest time increment, and a decrease in the period 
of the oscillatory response as the time increment 
is decreased. 

In figures 13 to 15, the time increment is kept 
constant, and the responses for idealization of the 
beam by 3 and 5 masses are plotted to bring out 
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Figure 11. Bending moment ratio at the center of the beam 
caused by impact at a tip, three-mass idealization. 
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Figure 12. Bending moment ratio at the center of the beam 
caused by impact at a tip, five-mass idealization. 
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Figure 13. Bending moment ratios for different mass ideali- 
zations with impact at tip, time increment At = ^mfi/12EI. 
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Figure 14. Bending moment ratios for different mass ideali- 
zations with impact at tip, time increment At= 1 /3^Jml z / 12EI . 
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Figure 15. Bending moment ratios for different mass ideali- 
zations with impact at tip, time increment At=l/5\mP/l£EI . 
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(b) Idealized beams 



Figure 16. Mass distribution considered for two beams normal 
to each other subjected to tip impact. 
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Figure 17. Bending moment at wing root for pair of beams 
representing airplane, each beam idealized by five masses. 
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■nensionless time, t/l2EI/ml 3 



Figure 18. Bending moment at fuselage root for pair of beams 
representing airplane, each beam idealized by Jive masses. 



the effect of increasing the number of masses. It 
can be seen that the period of the initial vibratory 
response decreases with an increase in the number 
of masses and that the damping increases. The 
maximum bending moment at the center is relatively 
unaffected. 

3.3. Pair of Beams Representing a Fuselage and 
Wing With Impact at Both Wing Tips 

An airplane can be approximated by two uniform 
beams at right angles to each other, one beam repre- 
senting the wing, the other the fuselage. Such a 
system, figure 16 (a), is investigated. Each beam 
is idealized as 5 masses, as shown in figure 16 (b). 
The suddenly applied constant loads, P, are consid- 
ered to act at each wing tip. The bending moments 
at the center of both the "fuselage" and "wing" beams 
are computed by using the same three time incre- 
ments in the numerical integration as were previ- 
ously used. The bending moment at the center of 
the struck beam (wing) is shown in figure 17 and 
that for the other beam (fuselage) is shown in figure 
18. The agreement between the bending moment 
ratios for the three values of time increment is only 
fair even for the two shortest time increments. As 
might have been expected from the analysis given 
of the single-degree-of -freedom system, the error 
due to using the longest time increment is evident 
in the marked damping of the oscillatory response. 

3.4. Discussion 

It is estimated that a time increment of about 1/30 
of the period in the fundamental mode is necessary 
for good accuracy in numerical-integration methods. 
This value is based on several factors. 

(1) The solution given in eq (26) corresponds to 
about 20 time increments per fundamental period. 
In this case the frequency coefficient, 0.961, is in 
error by about 4 percent. 

(2) The convergence indicated in the figures is 
not quite complete for the shortest time increment, 



which corresponds to 1/13.74 of the fundamental 
symmetrical mode period and only 1/5.00 of the 
lowest antisymmetrical mode period. A shorter 
time increment than this is certainly needed to give 
the response of the higher modes. 

(3) The recommended time increment is likely to 
be used for computing responses with slower con- 
vergence than the bending moment at the center of 
the beam and for structures with slower convergence 
than that of a uniform beam. 

(4) It is likely that computations by numerical 
methods will be carried out by automatic computing 
machines where the advantage of additional accuracy 
from a smaller time increment can be had with less 
penalty than when hand computing is used. 

It is noted from figures 8, 9, and 10 that, for each 
time increment considered, the response to an impact 
at the center is essentially the same for the 5- and 
7-mass idealizations of the beam but considerably 
higher for the .'3-mass system. For impacts at one 
end of the beam, figures 13, 14, and 15, the response 
shows less change in going from the 3- to the 5-mass 
idealizations than is the case for impact at the center. 
It would appear that a uniform beam should be 
approximated by a 5-mass idealization if errors in 
the bending moment at the center are to be small. 

In table 1, the maximum absolute values of 
bending moment ratios are given. For the beam 
with the impact at the center, the absolute value of 
the bending moment increased as the time inter- 
val for integration was decreased, but decreased 
as the number of masses used for the idealization of 
the distributed mass of the beam was increased. 
From this, it may well be that the errors caused by 
using a reasonable value of time interval for the 
numerical integration will be offset to some extent 
by the errors caused by using a finite number of 
masses. For the impact at a tip mass, the increase 
in bending moment with decrease in time interval is 
appreciable, but the decrease in bending moment with 
increase in number of masses is much less than for the 
case of impact at the center. 

Table 1. Maximum absolute value of bending moment at 
center of uniform beam struck with force P 



Time 


interval 


3-mass 


\MIPl\ 




5-mass 


7-mass 


Impact at center 






0. 355 
. 450 
.480 


0. 305 
. 383 
.410 


0. 295 
.368 
.400 


imPIUEI 


l/3Vm/3/l2S/_.__ 






l/5^ml*/12EI 


Impact at tip 






0. 355 
.450 
.480 


0. 333 
.448 
. 485 




^mWUEI 


l!Z\mlill2EI 






1/5V mPj\2EI 
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It is evident from figure 5 that the period of the 
response approaches the value for the exact solution 
for a beam idealized by 3 masses as the time interval 
in the integration is shortened. This same decrease 
in period with time increment is noted in figures 6 
and 7 for the 5- and 7- mass idealization of the beam. 
The period of the response also decreases as the num- 
ber of masses is increased, figures 7, 8, and 9. This 
latter result, at least qualitatively, verifies that ob- 
tained by Duncan [9]. He has obtained an exact 
solution showing that the periods of the computed 
modes of vibration of a uniform beam decrease as 
the number of masses chosen to represent the beam 
is increased, and that the magnitudes of the errors 
in the computed periods vary inversely as the square 
of the number of masses. 

4. Conclusions 

From the analyses of the responses of single-degree- 
of-freedom and many-degree-of-freedom systems to 
a suddenly applied force, it may be concluded that: 

a) Numerical integration methods give results 
that closely approximate the exact solution only if 
the increment of time between successive steps in 
the integration is small compared to all the natural 
periods of vibration of the system. In cases similar 
to those investigated where bending moment near 
the center is desired, a time increment of about %o 
of the fundamental period of the system will give 
good results. 

b) Of the methods investigated, the Houbolt 
method is the only one that gives convergent results 
for large time increments. In this method, errors 
introduced by large time increments result in the 
damping out of oscillatory response. 

c) The distributed mass of the beam can be con- 
sidered to be concentrated at relatively few mass 
points for computational purposes. In cases similar 
to those investigated, approximating a uniform beam 
by a 5-mass idealization gives good results. 

5. Appendix 

The detailed method of determining the response 
of a uniform beam to an impact load by numerical 
integration is given in this appendix for the case of 
a beam subjected to an impact at its tip. The case 
considered is that where the time increment is 
At—^lml 3 /l2EI. The responses for all the other 
cases presented in this report are obtained in a 
similar manner. 



Station 12 3 4 5 

Mass m/8 m/4 m/4 m/4 m/8 



Figure 19. Mass distribution for example considered in 
appendix. 



The mass of the uniform beam is considered to be 
distributed at five stations along its length, figure 19. 
The force, P, is applied at the tip, station 5. As 
used ia this report, an influence coefficient, 5 r , s is the 
displacement at station r due to a unit load at station 
s with the center of the beam, station 3, fixed. The 
influence coefficients are given in table 2. Denoting 
the displacement positive upward at station r by d r , 
the force at station r by F T , and the rotation at station 
3 by a, the displacements of the beam at stations 1, 
2, 4, and 5 are 



d l = dz + d n F 1 J r 8 l2 F 2 —al, 



d 2 = d 3 + h u Fi + d 22 F 2 — ol-1 



d A =d 3 + 544F4 + 8t 5 F 5 + a 



I 



(Al) 



d 5 =d3 J r 5 45 ^4 + h 5 F 5 + al ^ 



Substituting the values of the influence coefficients 
from table 2 into eq (Al), and solving for the forces 

Table 2. Influence coefficients 6 r , 8 for 5-mass idealized beam 
shoivn in figure 19 



\ 












\r 


1 


2 


3 


4 


5 


\ 












1 


W3EI 


5/3/48£7 











2 


bi*mEi 


J3/24EJ 











3 

















4 











PI2AEI 


5Z3/48£/ 


5 











5V/&EI 


WSE1 



in terms of the displacements at each station and of 
the rotation at station 3. 



F-2- 

F<= 

iv 



48£7, 



(-5(/, + 16c? 2 -llA+3a/), 

(-llds+Wdi-Bds-Sal), 

48Er/ , , , ,oi\ 
-jp- i 3da— 5d 4 +2dt+-^] 



' 7P 

48EI 

' 7P 

4SEI 



y (A2) 



From the condition that the sum of the forces equal 
zero, 

^1 + ^2 + ^3 + ^4+^5=0, (A3) 

and substituting from eq (A2), the value of F 3 is 
F l =^j^[M i -Ud i +\Qd i ~lUA 3iJ (A4) 



64 



From the requirement that the net moment equal _ rt _ d } , _ d 2 . - - d 4 A oe? rf 5 , AoN 

zerQj a=0.25y-1.5j+1.5-| — 0.25 y (A6) 



and using eq (A2) and (A4), 



5) 



Substituting the value of a from eq (A6) into eq (A2), 
and rewriting eq (A4): 



F 1 = ^|^(1.875rf 1 -4.25rf 2 -f^.00r/ : 3-0.75r/ 4 +0.125^ 5 ), 
^ 2 = ^-(-4.25r/ 1 + 11.50r/ 2 -11.00rf 3 + 4.50r/ 4 -0.75f/ 5 ), 
^3=1^^(3.00^- 11.00^2+ 16.00^3- 11.00^4+3.00^), 
F 4 = ^^ (-0.75^ + 4.50^-1 l.OOr/3+l .1.50f/ 4 -4.25</ 5 ), 



F*= 



7P 
48EI 



7P 



(0.125rfi-0.75rf 2 +3.00rf 3 — 4.25rf 4 +1.875rfs) 



(A7) 



The forces are the sum of the inertia loads due to 
the acceleration of the mass of the beam and the 
applied forces, 



F^-^l. 



F^-^l, 



F 3 =-jd 3 , 



F^-^K 



F b =-^d B +p 



where m is mass of uniform beam; d T , acceleration of 
beam at station r(r= 1,2, ... 5); and P, applied 
constant force at station 5. 

Using the difference eq (17), acceleration at station 
r can be expressed as 

2 5 4 1 



At 2 



At 2 



(A9) 



Equating the first of eq (A7) and (A8) and substi- 
tuting for (d T ) t from eq (A9), we obtain 



m 



— g^[2(d0|-8(il)l-Al + 4(rf0l-2A|-(rfl)l-8Af] = 



48£7 

7P 



[L875di— 4.25d 2 +3.00rf 3 — 0.75d 4 +0.125rfJ«. 

(A10) 



Similar relations are obtained by using the remain- 
ing eq (A 7) and (A8). 

Multiplying eq (AlO) and the corresponding equa- 
tions for the other stations through by SEl/PP to 
express the deflections as a dimensionless ratio, and 
transposing terms, the equations can be written in a 
form which, for station 1, reads 



0.125 (^p) =(«/,)«, (All) 



where 



4 (-w L,X w-),-uj (A12) 

Equation (All) and the equations similar to it can be 
expressed in matrix form as 

[A] ( ^ d j^ I ] ) ={J T }, (A13) 

where the matrix [A] is given in table 3. Solving eq 
(Al3) for the deflections gives, 

&r}r [A] ~ 1{J ' } " (A14) 
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Table 3. Matrix [A] of deflection coefficients for use in equation {A13') 



18-5 1 7ml3 


-4. 25 
7ra/3 

UJB +S«aF 

-11.00 

4.50 

-0.75 


3.00 
-11.00 
16-+- 7WZ3 
-11.00 

3.00 


-0.75 

4.50 

-11.00 

U5+ 96EIAt* 
-4.25 


0.125 

-.75 

3.00 

-4.25 

i q?R i 7w ' 3 
1 - 8/5+ 192£J/A«2 


L8 ' 5 ' 192EIAV 
-4.25 

3.00 

-0.75 

.125 



where [A] * is obtained as the inverse of [^4] with 
At=ymP/12EI. The matrix [^1] _1 is given in table 4. 
The first of this set of equations in conventional form 
is 

( 3 ^ / )=1.344937(J 1 ), + 0.589567(J 2 ),+ 

0.113294(J 3 );-0.113732(J 4 );-0.237485(J 5 )m (A 15) 
where, from eq (A 12), 

(J,),= 1.09375 Q^P) ^ -0.87500 Qyf 1 ) ^ 

+ 0.21875 (^l 7 ),^- (A16) 

Similar equations can be obtained for J at the other 
stations. 



Table 4. Inverted matrix [A]-* o/ taWe 3 iwtf M = ^Jm.l 3 /12EI 



1. 344937 


0. 589567 


0. 113294 


-0. 113732 


-0. 237485 


0. 589567 


. 519302 


. 305894 


. 079741 


-. 113732 


. 113294 


. 305894 


. 417773 


. 305894 


. 113294 


-. 113732 


. 079741 


. 305894 


. 519302 


. 589567 


-.237485 


-. 113732 


. 113294 


. 589567 


1.344937 



It is assumed that, prior to the application of the 
force, P, the displacements at all stations are zero. 
The computations are tabulated as shown in tables 
5 and 6. The J's are computed for time t, table 6, 
using eq (A16) for (Ji) t and similar equations for 
the other J's. The d's are computed for time t, 
table 5, using eq (A15) for (3diEI/Pl z ) t and similar 
equations for the other d's. From the d's at time t, 
the J's at time t+At are computed and, from these, 
the d's at time t-\-M. Thus a time history of the 
displacement ratios of the uniform beam due to the 
suddenly applied constant force is obtained. 

Since the step-by-step numerical procedure is 
dependent on previously computed values, it is 
desirable to have a check column which will indicate 
errors made in the computation. The deflections at 
each station were checked by computing the dis- 
placement, D, of the center of gravity of the beam 
from two different formulas. If these two values 
of D agreed, the computation of the deflections was 
assumed correct. The first equation for D t is: 



where 



m = X) m r 



(A 18) 



ZEI 
PP 



3EI 
PP 



^iS^-ff-J/ (A17) 



R_0.1U(?*«)+O.MO (?*«) + 
0.250 (^i) +0.250 ($*?) + 

0.125 (^ff); (A19) 



The other equation is obtained by substituting for 
m.EDIPP}^ [(3d 2 EI)/Pl% . . . [(3d 5 EI)/PP] t 
in eq (A 19), their equivalents in terms of the j's eq 
(A 14). This gave 

-pp- #,=0.28571375(Ji), + 0.28571362(J 2 ),+ 

0. 28571 375(J 3 )* + 0. 28571 362(J 4 )* + 
0.28571375(J 5 ),. (A20) 

The small differences in the coefficients in eq (A20) 
result from rounding-off error. For a good error 
indication, the coefficients, as given, are needed. 
The values of D from eq (A19) and (A20) are given 
in columns (7) and (8) of table 5. 

The J's are checked by comparing their sums as 
determined from two different formulas. The first 
formula is 

(S Jr)=(Jl)t + (J2)t + Wt + Wt + (J^ (A21) 

The second formula is obtained by substituting in eq 
(A21) the value of (Ji) t in eq (A16) and correspond- 
ing values for the other J's. Then, making use of eq 

(A19), 

(S Jr)=8.75(Z?),_ A «-7.00(Z>),_ 2A ,+ 

1.75(Z>),_ 3 A*+0.4375. (A22) 

The values computed from eq (A21) and (A22) are 
given in columns (7) and (8), respectively, of table 
6. The good agreement between the check columns 
in tables 5 and 6 indicate that the computation is 
free of numerical error. 
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Table 5. Displacement ratios at stations along a uniform beam subjected to an impact at its tip 



• 


2 


3 


4 


5 


6 


7 


8 


9 


Time 
ratio 

- I12EI 




Displacement ratio S(l T EI/Pl s 


at stations— 




SEID/PP 


Bending 

moment 

ratio at 

center, M/Pl 


1 


2 


3 


4 


5 


Eq (A19) 


Eq (A 20) 


-3 

-2 

-1 



1 






-0. 103900 
-. 525637 






-0. 049758 
-. 168315 






0. 049566 
. 294203 






0. 257936 
. 967785 








0.588410 

1.83828 






0. 125000 
. 437499 






0. 125000 
. 437499 







. 153654 
. 309800 


2 
3 

4 
5 
6 


- 1. 40216 
-2. 72217 
-4. 42824 
-6. 50884 
-8. 99357 


-.348058 
-. 596199 
-. 927607 
-1.34718 
-1.85109 


. 820339 
1.61553 
2. 64076 
3. 88834 
5. 37772 


2. 20772 
4.00162 
6. 35719 
9. 28142 
12. 7830 


3. 79212 
6. 55516 
10. 2248 
14. 8320 
20. 3579 


. 968745 
1. 73436 
2. 74216 
3. 99603 
5. 49794 


. 908745 
1.73436 

2. 74216 

3. 99603 
5. 49794 


. 322636 
. 253868 
. 213325 
. 227965 
. 257022 


7 
8 
9 
10 
11 


-11.9070 
-15.2485 
-19.0059 
-23. 1735 
-27. 7547 


-2.43581 
-3. 10178 
-3. 85121 
-4. 68514 
-5. 60308 


7. 12514 
9. 13013 
11.3845 
13. 8843 
16. 6319 


16. 8678 
21. 5368 
26. 7885 
32. 6223 
39. 0387 


26. 7835 
34.1121 
42. 3569 
51.5247 
61. 6123 


7. 24885 
9. 24923 
11.4993 
13. 9993 
16. 7491 


7. 24885 
9. 24923 
11.4993 
13. 9993 
16. 7491 


. 265110 
. 254424 
. 244509 
. 244722 
. 249966 


12 
13 
14 
15 
16 


-32. 7545 
-38. 1739 
-44.0108 
-50. 2639 
-56.9335 


-6. 60427 
-7.68860 
-8. 85645 
-10.1081 
-11.4436 


19.6304 
22. 8806 
26. 3810 
30. 1306 
34. 1295 


46. 0388 
53. 6226 
61. 7899 
70. 5404 
79. 8742 


72. 6149 
84. 5319 
97. 3652 
111.117 
125. 786 


19. 7488 
22. 9984 
26. 4979 
30. 2473 
34. 2465 


19. 748S 
22. 9984 
26. 4979 
30. 2473 
34. 2465 


. 252659 
. 251377 
. 249285 
. 248816 
. 249636 


17 
18 
19 

20 


-64.0204 
-71.5252 
-79.4478 

-87. 7877 


-12.8628 
-14.3057 
-15.9524 
- 1 7. 6230 


38. 3783 
42. 8773 
47. 6262 
52. 6249 


89. 7915 
100. 292 
111.377 
123.044 


141.372 

157.874 
175.294 
193.631 


38. 4957 
42. 9946 
47. 7434 
52. 7420 


38. 4957 
42. 9946 
47. 7434 
52. 7420 


. 250307 
. 25022S 
. 249824 
. 249634 



Table 6. Values of J at stations on beam 



1 


2 


3 


4 


5 





7 


8 


Time 
ratio, 

, \\2EI 




Values 


of J at station 


S— 




eq (A 22) 


1 


2 


3 


4 


5 


-3 
-2 

-1 


1 







-0. 113640 






-0. 108845 






0. 108426 






0. 564234 








0. 437500 

1.08107 






0. 437500 

1. 53125 






0. 437500 

1. 53125 


2 
3 
4 
5 
6 


-. 484003 
-1.09641 
-1.85547 
-2. 76821 
-3. 83981 


-.281114 
-.488595 
-. 768722 
-1. 13807 
-1.58449 


. 556827 
1.30132 
2. 22709 
3. 30839 
4.59119 


1. 66564 
3. 24861 
5. 31345 
7. 86939 
10. 9287 


1. 93326 
3. 10535 
4. 69122 
6. 71464 
9. 14720 


3.39001 
6. 07028 
9. 59757 
13. 9861 
19. 2428 


3.39061 
fi. 07028 
9. 59757 
13. 9861 
19.2428 


7 
8 
9 
10 
11 


-5.11017 
-6. 57770 
-8 22675 
- 10. 0499 
-12.0515 


-2. 09753 
-2. 67831 
-3. 33234 
-4. 06206 
-4. 86616 


6.11451 
7. 87638 
9. 85591 
12.0431 
14. 4435 


14. 5016 
18. 5887 
23. 1857 
28. 2902 
33. 9036 


11.9626 
15. 1633 
18. 7653 
22. 7763 
27. 1923 


25.3710 
32. 3724 
40. 2478 
48. 9975 
58. 6219 


25. 3710 
32. 3724 
40. 2478 
48. 9975 
58. 6219 


12 
13 
14 
15 
16 


-14.2374 
-16.6091 
-19. 1638 
-21. 8998 
-24. 8172 


-5. 74265 
-6. 69120 
-7. 71269 
-8. 80780 
-9. 97647 


17. 0654 
19.9102 
22. 9745 
26. 2557 
29. 7542 


40. 0282 
46. 6643 
53. 8110 
61. 4678 
69. 6347 


32. 0074 
37. 2204 
42. 8339 

48. 8498 
55. 2680 


69. 1209 
80. 4946 
92. 7428 
105. 866 
119. 863 


69. 1209 
80. 4940 
92. 7428 
105. 866 
119.863 


17 
18 
19 
20 


-27.9174 
-31. 2008 
-34. 6671 
-38. 3159 


- 11. 2184 
-12.5334 
-13.9217 
-15.3834 


33.4714 
37. 4081 
41. 5637 
45. 9376 


78. 3122 
87. 5005 
97. 1993 
107. 409 


62. 0872 
69. 3071 
76. 9279 
84. 9502 


134. 735 
150. 481 
167. 102 
184. 597 


134. 735 

150. 481 
167. 102 
184. 597 
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The forces are determined from the displacements 
of the masses at the five stations on the beam by use 
of eq (A7). The forces, multiplied by thier distances 
from the center of the beam determine the bending 
moment at the center of the beam. Designating the 
bending moment at the center by M, 



M=F,l+F 2 l - 
^M^(— 0.25rfi+1.50d 2 -2.50d 8 + 



(A23) 



ll 2 



1.50d 4 — 0.25d fi ) (A24) 



g~0.5714286 5*^+8.4285714 ^ / - 

5.7142857 3 ^ / + 3.4285714 ^~- 

0.5714286 ^§~ (A25) 

Values of the bending moment ratio computed 
from eq (A25) are given in column 9, table 5. They 
are plotted in figure 12 for At=^Jml 3 /12EI. 



The authors thank Airs. L. W. Roberson for com- 
puting the responses of the uniform beams to impact 
loads. 
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